Abstract
Part of theR for Artists and Designers course at the School of Foundation Studies, Srishti Manipal Institute of Art, Design, and Technology, Bangalore.
At the end of this Lab session, we should: - know the types and structures of spatial data and be able to work with them - understand the basics of modern spatial packages in R - be able to specify and download spatial data from the web, using R - plot static and interactive maps using ggplot, tmap and leaflet packages - add symbols and markers for places and regions of our own interest in these maps. - see directions for further work (e.g. maps + networks together)
rgdal is R’s interface to the “Geospatial Abstraction Library (GDAL)” which is used by other open source GIS packages such as QGIS and enables R to handle a broader range of spatial data formats.
We will take small steps in making maps using just two of the several map making packages in R.
The steps we will use are:
prettymapr or similar..)osmdataosmplot, tmap and also with leaflet.leaflet using a variety of map data providers. Note: tmap can also do interactive maps which we will explore also.Bas. Onwards and Map-wards!!
Let’s get BLR data into R and see if we can plot an area of interest. Then we can order on Swiggy…never mind.
Where is my home? Specify a “bounding box” first, using a rough longitude latitude info, or using a place name and searching for the long/lat:
# BLR Bounding Box
bbox <- osmplotr::get_bbox(c(77.56,12.93,77.63,12.96))
bbox_l <- osmdata::getbb("Bangalore, India")
bbox_p <- prettymapr::searchbbox("Bangalore")
## Using default API key for pickpoint.io. If batch geocoding, please get your own (free) API key at https://app.pickpoint.io/sign-up
bbox
## min max
## x 77.56 77.63
## y 12.93 12.96
bbox_l
## min max
## x 77.46010 77.78405
## y 12.83401 13.14366
bbox_p # identical with bbox_l
## min max
## x 77.46010 77.78405
## y 12.83401 13.14366
# Get Map data
dat_B <- extract_osm_objects (key = "building", bbox = bbox)
## Issuing query to Overpass API ...
## Rate limit: 0
## Query complete!
## converting OSM data to sf format
dat_H <- extract_osm_objects (key = 'highway', bbox = bbox)
## Issuing query to Overpass API ...
## Rate limit: 0
## Query complete!
## converting OSM data to sf format
dat_P <- extract_osm_objects (key = 'park', bbox = bbox)
## Issuing query to Overpass API ...
## Rate limit: 0
## Query complete!
## converting OSM data to sf format
dat_G <- extract_osm_objects (key = 'landuse', value = 'grass', bbox = bbox)
## Issuing query to Overpass API ...
## Rate limit: 0
## Query complete!
## converting OSM data to sf format
dat_T <- extract_osm_objects (key = 'natural', value = 'tree', bbox = bbox)
## Issuing query to Overpass API ...
## Rate limit: 0
## Query complete!
## converting OSM data to sf format
#Useful keys include building, waterway, natural, grass, park, amenity, shop, boundary, and highway.
blr_map <- osm_basemap(bbox = bbox, bg = "gray20") %>%
add_osm_objects(., dat_B, col = "gray40") %>%
add_osm_objects(., dat_H, col = "gray80") %>%
add_osm_objects(., dat_G, col = "darkseagreen1") %>%
add_osm_objects(., dat_P, col = "darkseagreen") %>%
add_osm_objects(., dat_T, col = "green")
print_osm_map(blr_map)
#print_osm_map (map, filename = "Blr.jpeg", device = "jpeg", width = 10, units = "cm")
#print_osm_map (map, filename = "map.png", width = 10, units = "in", dpi = 300)
blr_map_2 <- osm_basemap(bbox = bbox, bg = "gray20") %>%
add_osm_objects(., dat_B, col = "gray40") %>%
add_osm_objects(., dat_H, col = "gray80") %>%
add_osm_objects(., dat_G, col = "darkseagreen1") %>%
add_osm_objects(., dat_P, col = "darkseagreen") %>%
add_osm_objects(., dat_T, col = "green")
print_osm_map(blr_map)
class(dat_B)
## [1] "sf" "data.frame"
names(dat_B)
## [1] "osm_id" "name"
## [3] "Friary" "addr.city"
## [5] "addr.district" "addr.full"
## [7] "addr.housename" "addr.housenumber"
## [9] "addr.postcode" "addr.state"
## [11] "addr.street" "addr.suburb"
## [13] "air_conditioning" "alt_name"
## [15] "amenity" "area"
## [17] "atm" "barrier"
## [19] "brand" "brand.wikidata"
## [21] "brand.wikipedia" "building"
## [23] "building.colour" "building.levels"
## [25] "building.levels.underground" "bus"
## [27] "capacity" "clothes"
## [29] "club" "construction"
## [31] "cuisine" "delivery"
## [33] "denomination" "description"
## [35] "designation" "ele"
## [37] "email" "emergency"
## [39] "government" "healthcare"
## [41] "healthcare.speciality" "height"
## [43] "internet_access" "landmark"
## [45] "landuse" "layer"
## [47] "leisure" "level"
## [49] "lit" "man_made"
## [51] "max_age" "min_age"
## [53] "name.en" "name.kn"
## [55] "note" "office"
## [57] "official_name" "opening_hours"
## [59] "operator" "operator.type"
## [61] "organic" "outdoor_seating"
## [63] "payment.cash" "payment.debit_cards"
## [65] "phone" "public_transport"
## [67] "religion" "screen"
## [69] "service" "service_times"
## [71] "shop" "smoking"
## [73] "source" "takeaway"
## [75] "tourism" "website"
## [77] "wheelchair" "wikidata"
## [79] "wikipedia" "geometry"
class(dat_B$geometry)
## [1] "sfc_POLYGON" "sfc"
nrow(dat_B)
## [1] 39699
We can create areas of interest around the map. Say based on metro stations or your restaurants etc.
my_area <-
bbox %>%
zoombbox(factor = 16, offset = c(0,0))
my_area
## min max
## x 77.59281 77.59719
## y 12.94406 12.94594
bbox
## min max
## x 77.56 77.63
## y 12.93 12.96
my_area_in_blr <-
sp::SpatialPoints(coords = cbind(
c(my_area["x", "min"], my_area["x", "max"],
my_area["x", "max"], my_area["x", "min"]),
c(my_area["y", "min"], my_area["y", "min"],
my_area["y", "max"], my_area["y", "max"])
)
)
my_area_in_blr
## class : SpatialPoints
## features : 4
## extent : 77.59281, 77.59719, 12.94406, 12.94594 (xmin, xmax, ymin, ymax)
## crs : NA
blr_map %>% add_osm_groups(
.,
dat_B,
groups = my_area_in_blr,
cols = "yellow",
bg = "gray40",
colmat = FALSE
) %>%
add_osm_objects(., dat_B, col = "red", size = 0.2) %>%
print_osm_map(.)
rnaturalearth and tmapdata("World")
data("metro")
head(metro, n = 3)
## Simple feature collection with 3 features and 12 fields
## Geometry type: POINT
## Dimension: XY
## Bounding box: xmin: 3.04197 ymin: -8.83682 xmax: 69.17246 ymax: 36.7525
## Geodetic CRS: WGS 84
## name name_long iso_a3 pop1950 pop1960 pop1970 pop1980 pop1990
## 2 Kabul Kabul AFG 170784 285352 471891 977824 1549320
## 8 Algiers El Djazair (Algiers) DZA 516450 871636 1281127 1621442 1797068
## 13 Luanda Luanda AGO 138413 219427 459225 771349 1390240
## pop2000 pop2010 pop2020 pop2030 geometry
## 2 2401109 3722320 5721697 8279607 POINT (69.17246 34.52889)
## 8 2140577 2432023 2835218 3404575 POINT (3.04197 36.7525)
## 13 2591388 4508434 6836849 10428756 POINT (13.23432 -8.83682)
tmap_mode("plot")
## tmap mode set to plotting
tm_shape(World) +
tm_polygons("HPI") +
tm_shape(metro) +
tm_bubbles(size = "pop2020", col = "red")
tmap_mode("view")
## tmap mode set to interactive viewing
tm_shape(World) +
tm_polygons("HPI") +
tm_shape(metro) +
tm_bubbles(size = "pop2020", col = "red") +
tm_tiles("Stamen.TonerLabels")
## Legend for symbol sizes not available in view mode.
india <-
ne_states(country = "india",
geounit = "india",
returnclass = "sf")
india_neighbours <-
ne_states(country= (c("india", "sri lanka", "pakistan", "afghanistan", "nepal","bangladesh")))
names(india)
## [1] "featurecla" "scalerank" "adm1_code" "diss_me" "iso_3166_2"
## [6] "wikipedia" "iso_a2" "adm0_sr" "name" "name_alt"
## [11] "name_local" "type" "type_en" "code_local" "code_hasc"
## [16] "note" "hasc_maybe" "region" "region_cod" "provnum_ne"
## [21] "gadm_level" "check_me" "datarank" "abbrev" "postal"
## [26] "area_sqkm" "sameascity" "labelrank" "name_len" "mapcolor9"
## [31] "mapcolor13" "fips" "fips_alt" "woe_id" "woe_label"
## [36] "woe_name" "latitude" "longitude" "sov_a3" "adm0_a3"
## [41] "adm0_label" "admin" "geonunit" "gu_a3" "gn_id"
## [46] "gn_name" "gns_id" "gns_name" "gn_level" "gn_region"
## [51] "gn_a1_code" "region_sub" "sub_code" "gns_level" "gns_lang"
## [56] "gns_adm1" "gns_region" "min_label" "max_label" "min_zoom"
## [61] "wikidataid" "name_ar" "name_bn" "name_de" "name_en"
## [66] "name_es" "name_fr" "name_el" "name_hi" "name_hu"
## [71] "name_id" "name_it" "name_ja" "name_ko" "name_nl"
## [76] "name_pl" "name_pt" "name_ru" "name_sv" "name_tr"
## [81] "name_vi" "name_zh" "ne_id" "geometry"
names(india_neighbours)
## [1] "featurecla" "scalerank" "adm1_code" "diss_me" "iso_3166_2"
## [6] "wikipedia" "iso_a2" "adm0_sr" "name" "name_alt"
## [11] "name_local" "type" "type_en" "code_local" "code_hasc"
## [16] "note" "hasc_maybe" "region" "region_cod" "provnum_ne"
## [21] "gadm_level" "check_me" "datarank" "abbrev" "postal"
## [26] "area_sqkm" "sameascity" "labelrank" "name_len" "mapcolor9"
## [31] "mapcolor13" "fips" "fips_alt" "woe_id" "woe_label"
## [36] "woe_name" "latitude" "longitude" "sov_a3" "adm0_a3"
## [41] "adm0_label" "admin" "geonunit" "gu_a3" "gn_id"
## [46] "gn_name" "gns_id" "gns_name" "gn_level" "gn_region"
## [51] "gn_a1_code" "region_sub" "sub_code" "gns_level" "gns_lang"
## [56] "gns_adm1" "gns_region" "min_label" "max_label" "min_zoom"
## [61] "wikidataid" "name_ar" "name_bn" "name_de" "name_en"
## [66] "name_es" "name_fr" "name_el" "name_hi" "name_hu"
## [71] "name_id" "name_it" "name_ja" "name_ko" "name_nl"
## [76] "name_pl" "name_pt" "name_ru" "name_sv" "name_tr"
## [81] "name_vi" "name_zh" "ne_id"
#ind <- metro$iso_a3 == "IND"
#metro_ind1 <- metro[ind,]
metro_ind2 <- subset(metro, iso_a3 == "IND")
metro_neighbours <- metro %>% dplyr::filter(iso_a3 %in% c("IND","PAK", "LKA", "BGD","NPL"))
tm_shape(World %>% dplyr::filter(iso_a3 %in% c("IND", "AFG", "PAK", "NPL", "BGD", "LKA"))) +
tm_borders() +
tm_shape(india) +
tm_polygons("name") +
tm_shape(metro_ind2)+
tm_dots(size = "pop2020") +
tm_layout(legend.outside = TRUE,legend.outside.position = "right") +
tm_credits("Geographical Boundaries are not accurate",size = 0.5,position = "right") +
tm_compass(position = c("right", "top")) +
tm_scale_bar(position = "left") +
tmap_style("watercolor") #cobalt #gray #white #col_blind #beaver #classic #watercolor #albatross #bw
## tmap style set to "watercolor"
## other available styles are: "white", "gray", "natural", "cobalt", "col_blind", "albatross", "beaver", "bw", "classic"
## Credits not supported in view mode.
## Compass not supported in view mode.
## Warning: Number of levels of the variable "name" is 35, which is
## larger than max.categories (which is 30), so levels are combined. Set
## tmap_options(max.categories = 35) in the layer function to show all levels.
## Legend for symbol sizes not available in view mode.
tm_shape(india_neighbours) +
tm_polygons("name") +
tm_shape(metro_neighbours) +
tm_dots(size = "pop2020") +
tm_layout(legend.outside = TRUE,legend.outside.position = "right") +
tmap_options(max.categories = 10) +
tm_credits("Geographical Boundaries are not accurate",size = 0.5,position = "center")
## Credits not supported in view mode.
## Warning: Number of levels of the variable "name" is 112, which is
## larger than max.categories (which is 10), so levels are combined. Set
## tmap_options(max.categories = 112) in the layer function to show all levels.
## Legend for symbol sizes not available in view mode.
tmap_mode("plot")
## tmap mode set to plotting
tm_basemap("Stamen.Watercolor") +
tm_shape(metro, bbox = "India") +
tm_dots(col = "red",
# user-chosen group name for layers
group = "Metropolitan Areas") +
tm_shape(World) +
tm_borders() +
tm_tiles(server = "Stamen.TonerLabels",group = "Labels") + # ADDS LABELS!!!
tm_graticules()